%% read data
fname = 'D:\Manoj\tsunami\obs_data\paf20041224d.min';
data_struct = intermagnet(fname , 'HDZF');
paf_z = data_struct.DATA(:,3);
paf_h = data_struct.DATA(:,1);
paf_fday = data_struct.FDAY;

fname = 'D:\Manoj\tsunami\obs_data\paf20041225d.min';
data_struct = intermagnet(fname , 'HDZF');
paf_z(end+1:end+1440) = data_struct.DATA(:,3);
paf_h(end+1:end+1440) = data_struct.DATA(:,1);
paf_fday(end+1:end+1440) = data_struct.FDAY;

fname = 'D:\Manoj\tsunami\obs_data\paf20041226d.min';
data_struct = intermagnet(fname , 'HDZF');
paf_z(end+1:end+1440) = data_struct.DATA(:,3);
paf_h(end+1:end+1440) = data_struct.DATA(:,1);
paf_fday(end+1:end+1440) = data_struct.FDAY;


fname = 'D:\Manoj\tsunami\obs_data\ams20041224d.min';
data_struct = intermagnet(fname , 'HDZF');
ams_z = data_struct.DATA(:,3);
ams_h = data_struct.DATA(:,1);
ams_fday = data_struct.FDAY;

fname = 'D:\Manoj\tsunami\obs_data\ams20041225d.min';
data_struct = intermagnet(fname , 'HDZF');
ams_z(end+1:end+1440) = data_struct.DATA(:,3);
ams_h(end+1:end+1440) = data_struct.DATA(:,1);
ams_fday(end+1:end+1440) = data_struct.FDAY;

fname = 'D:\Manoj\tsunami\obs_data\ams20041226d.min';
data_struct = intermagnet(fname , 'HDZF');
ams_z(end+1:end+1440) = data_struct.DATA(:,3);
ams_h(end+1:end+1440) = data_struct.DATA(:,1);
ams_fday(end+1:end+1440) = data_struct.FDAY;


fname = 'D:\Manoj\tsunami\obs_data\czt20041224d.min';
data_struct = intermagnet(fname , 'HDZF');
czt_z = data_struct.DATA(:,3);
czt_h = data_struct.DATA(:,1);
czt_fday = data_struct.FDAY;

fname = 'D:\Manoj\tsunami\obs_data\czt20041225d.min';
data_struct = intermagnet(fname , 'HDZF');
czt_z(end+1:end+1440) = data_struct.DATA(:,3);
czt_h(end+1:end+1440) = data_struct.DATA(:,1);
czt_fday(end+1:end+1440) = data_struct.FDAY;

fname = 'D:\Manoj\tsunami\obs_data\czt20041226d.min';
data_struct = intermagnet(fname , 'HDZF');
czt_z(end+1:end+1440) = data_struct.DATA(:,3);
czt_h(end+1:end+1440) = data_struct.DATA(:,1);
czt_fday(end+1:end+1440) = data_struct.FDAY;


%% wavelet analysis

load d:\manoj\tsunami\obs_data\paf_czt_ams_dec24_26_2004
b1 = min(ams_fday):(1/24):max(ams_fday);
sp=spline(b1,ams_h'/spline(b1,eye(length(b1)),ams_fday'));
ams_h_f=ams_h-ppval(ams_fday,sp);

sp=spline(b1,ams_z'/spline(b1,eye(length(b1)),ams_fday'));
ams_z_f=ams_z-ppval(ams_fday,sp);

sp=spline(b1,paf_h'/spline(b1,eye(length(b1)),paf_fday'));
paf_h_f=paf_h-ppval(paf_fday,sp);

sp=spline(b1,paf_z'/spline(b1,eye(length(b1)),paf_fday'));
paf_z_f=paf_z-ppval(ams_fday,sp);

sp=spline(b1,czt_h'/spline(b1,eye(length(b1)),czt_fday'));
czt_h_f=czt_h-ppval(czt_fday,sp);

sp=spline(b1,czt_z'/spline(b1,eye(length(b1)),czt_fday'));
czt_z_f=czt_z-ppval(czt_fday,sp);


%clear ams_z ams_h paf_z paf_h;
%Sc = 1:5:100;
%Sc = 2:5:200;
Sc = 1:1:100;
%waven = 'morl';
waven = 'cgau4';
perd = 1./scal2frq(Sc,waven,1);


pafzw = cwt(paf_z_f,Sc,waven);%Complex wavelet transfor
%clear paf_z_f;
amszw = cwt(ams_z_f,Sc,waven);%Complex wavelet transform
%clear ams_z_f;
amshw = cwt(ams_h_f,Sc,waven);%Complex wavelet transform
%clear ams_h_f;
pafhw = cwt(paf_h_f,Sc,waven);%Complex wavelet transformf
%clear paf_h_f;
cztzw = cwt(czt_z_f,Sc,waven);%Complex wavelet transformf
czthw = cwt(czt_h_f,Sc,waven);%Complex wavelet transformf

Hxy = amshw .* conj(pafhw);
L = abs(Hxy) > 1000;
Hxy(L) = 1000;
ss = sqrt(abs(Hxy));
cc = max(max(ss)) - abs(ss);
weight = (cc./max(max(cc)));
%clear cc ss L amshw pafhw paf_fday;


%% plot z
Eq = datenum(2004,12,26,0,58,0);
Ams_arrival = datenum(2004,12,26,7,58,0);
Paf_arrival = datenum(2004,12,26,9,58,0);
subplot(311);
imagesc(ams_fday,perd,abs(amszw));
hold on;
datetick('x','keeplimits');
line([Ams_arrival Ams_arrival],[100 0],'color',[0.7 0.7 0.7],'LineWidth',2)
caxis([0,20]);
subplot(312);
imagesc(ams_fday,perd,abs(pafzw));
hold on;
line([Paf_arrival Paf_arrival],[100 0],'color',[0.7 0.7 0.7],'LineWidth',2)
datetick('x','keeplimits');
caxis([0,20]);
subplot(313);
imagesc(ams_fday,perd,abs(cztzw));
datetick('x','keeplimits');
caxis([0,20]);

%% plot h
Eq = datenum(2004,12,26,0,58,0);
Ams_arrival = datenum(2004,12,26,7,58,0);
Paf_arrival = datenum(2004,12,26,9,58,0);
subplot(311);
imagesc(ams_fday,perd,abs(amshw));
hold on;
datetick('x','keeplimits');
line([Ams_arrival Ams_arrival],[100 0],'color',[0.7 0.7 0.7],'LineWidth',2)
caxis([0,20]);
subplot(312);
imagesc(ams_fday,perd,abs(pafhw));
hold on;
line([Paf_arrival Paf_arrival],[100 0],'color',[0.7 0.7 0.7],'LineWidth',2)
datetick('x','keeplimits');
caxis([0,20]);
subplot(313);
imagesc(ams_fday,perd,abs(czthw));
datetick('x','keeplimits');
caxis([0,20]);
